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The parallel-tempering method has been applied to numerically study the thermo- 
dynamic behavior of a three-dimensional disordered antiferromagnetic Ising model 
with random fields at spin concentrations corresponding to regions of both weak and 
strong structural disorder. An analysis of the low-temperature behavior of the model 
convincingly shows that in the case of a weakly disordered samples there is realized 
an antiferromagnetic ordered state, while in the region of strong structural disorder 
the effects of random magnetic fields lead to the realization of a new phase state 
of the system with a complex domain structure consisting of antiferromagnetic and 
ferromagnetic domains separated by regions of a spin-glass phase and characterized 
by a spinglass ground state. 

I. INTRODUCTION 

Studying the critical behavior of disordered systems with quenched structural defects 
is of a great theoretical and experimental interest since the majority of real solids contain 
quenched structural defects whose presence influences their thermodynamic characteristics 
and, in particular, substantially affects the behavior of the systems during phase transitions. 
It is known that the effect of quenched impurities manifests itself in the form of random 
disturbances of the local temperature in ferro and antiferromagnetic systems in the absence 
of an external magnetic field or in the form of random magnetic fields in antiferromagnetic 
systems in a uniform magnetic field. 

The investigations performed showed that in the first case the presence of quenched 
nonmagnetic impurity atoms changes the properties of only anisotropic Ising like magnets 
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in phase transitions in which the heat capacity in the purely uniform state diverges at a 
critical temperature. In the opposite case of magnets described by the XY model or by the 
isotropic Heisenberg model, the presence of an impurity does not affect their behavior during 
phase transitions. The critical properties of the disordered Ising model has recently been 
studied in numerous works (see references in For the dilute Ising-like systems, good 
agreement was obtained between the theoretical calculations and the results of experiments 
and computer simulation by the Monte Carlo method. An opposite situation is observed 
for the magnetic systems with a disorder of the random magnetic field type. In spite of 
numerous studies, which began since 1975 when this type of disorder was described for the 

n 

first time [2|], only scarce reliable information on the behavior of these systems exists at 
present jsj. In particular, the nature of the phase transition in the three-dimensional Ising 
model with random fields still remains unexplained, and the theoretical results obtained 
contradict the experimental data. The comparison of theoretical predictions with the results 
of experimental studies is hindered by the difficulties of achieving the equilibrium state in 
such systems because of their anomalously slow relaxation properties. According to some 
data, this is a first-order transition up to very low values of random fields; according 



to other, this is a second-order transition 



For describing the influence of random fields on the behavior of magnetic systems, two 
qualitatively equivalent models, namely, the ferromagnetic random-field Ising model (RFIM) 



9] and the disordered antiferromagnetic Ising model in an external uniform field (DAFF) 



1 01 ] are used. The real magnetic systems with the effects of random fields are antiferro- 
magnets with quenched impurities of nonmagnetic atoms; in the behavior of these systems, 
effects of a ferromagnetic interaction of next-nearest atoms manifest themselves along with 
the antiferromagnetic interaction of nearestneighbor atoms. The structure of an antiferro- 
magnet can be represented as several ferromagnetic sublattices inserted into each other in 
such a way that the total magnetization of the antiferromagnet remains equal to zero despite 
the fact that at a temperature lower than the Neel point there occurs a magnetic ordering 
within each ferromagnetic sublattice. As examples of two-sublattice antiferromagnets, NiO, 
MnO, Fe2C>3, MnF 2 , and some other can be mentioned. As the examples of the realization 
of disordered systems with random magnetic fields, uniaxial Ising-like antiferromagnets such 
as MnF 2 and FeF 2 with impurities of Zn atoms in an external magnetic field can be taken 
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In [12l . Il3j |. we for the first time showed, using Monte Carlo computer simulation of the 
thermodynamic behavior of the disordered antiferromagnetic Ising model with effects of ran- 
dom magnetic fields, that, in weakly disordered systems with a spin concentration higher 
than the threshold for impurity percolation, a second-order phase transition from the para- 
magnetic into the antiferromagnetic state is realized. In strongly disordered systems with a 
spin concentration lower than this threshold value, a first-order transition from the param- 
agnetic into a mixed state, which is characterized by a complex domain structure consisting 
of antiferromagnetic and ferromagnetic domains separated by regions of a spinglass phase, 
occurs in the system. It was shown that, with a reduction in the spin concentrations and 
an increase in the magnitude of the external magnetic field, a decrease in the number and 
sizes of antiferromagnetic domains occurs in the system, and an increase in the number and 
sizes of ferromagnetic domains along with a reduction of the relative volume of the spin-glass 
phase is observed. It is shown that in this region of spin concentrations the effects of random 
magnetic fields lead to a change from an antiferromagnetic ground state to a spin-glass state. 

In this work, we consider the same antiferromagnetic Ising model (with a spin concentra- 
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tion p = 0.5 corresponding to the region of strong structural disordering) as in 12. 131. The 
aim is to obtain an additional confirmation of the existence in such systems of a spin-glass 
ground state and of a complex domain structure by the realization and application, for its 
numerical study, of the algorithm of the parallel-tempering method which was developed 
specially for studying the thermodynamics of spin glasses. To compare the results obtained, 
we also investigated the model at p = 0.9 corresponding to the region of weak disorder. 



The disordered two-sublattice antiferromagnetic Ising model was defined as a system of 
spins with a concentration p connected with N = pL 3 sites of a cubic lattice with periodic 
boundary conditions. The Hamiltonian of the model under consideration has the form 



where Oi = ±1, /i is the Bohr magneton, J\ = 1 and J<i = —1/2 characterize the antifer- 
romagnetic interaction of spins with nearest neighbors and their ferromagnetic interaction 
with next-nearest neighbors, respectively; and h is the intensity of a uniform magnetic field. 



II. MODEL 




(1) 



The random variables Pi and pj are described by the distribution function 



P(Pi)=p5( Pl -l) + (l-p)5(p l ) 



(2) 



and characterize the quenched nonmagnetic impurity atoms distributed over the lattice sites 
(empty sites). 

The simulation of statistical properties of this model made it possible to determine ther- 
modynamic values such as the total magnetization 



the staggered magnetization M stg = Mi — M2 (Mi and M2 - are the magnetizations of the 
sublattices), and the spin-glass order parameter 



where the exponents a and (3 characterize different replicas of the disordered system that 
are simulated simultaneously at one and the same temperature and have different initial 
configurations. In expressions ([3]) and (HI), the angular brackets designate statistical averag- 
ing realized for each impurity configuration of the system, and the brackets mean averaging 
over different impurity configurations. 

The quantities M, M stg , and q at p, characterize different types of magnetic ordering of a 
strongly disordered system that can appear in it in the low-temperature phase. Along with 
these magnetic thermodynamic quantities, measurements of the heat capacity as a thermal 
characteristic of occurring phase transformations in the system were realized. 



As is known, the spin-glass state is characterized by the presence of a large number 
of metastable energy states separated by potential barriers. The number of metastable 
states exponentially grows with increasing number of spins, which strongly complicates 
the numerical simulation of such systems. In spin glasses, there is a problem of reach- 
ing an equilibrium state because of the existence of high energy barriers which separate 
local energy minima. At sufficiently low temperatures, the system can never leave a lo- 
cal energy minimum even if the corresponding state is globally unstable. This feature 




(3) 





III. PARALLEL-TEMPERING METHOD 



makes impossible obtaining physical characteristics for the magnetic materials which con- 
tain, as in our case, a spin-glass phase when using standard Monte Carlo algorithms. 
There arises a need for an improvement or a modi- 
fication of the algorithms used for simulation. One 
of the algorithms which made it possible to solve 
this problem became the parallel-temperature al- 



gorithm 



The parallel-tempering method is an extension 
of the usual Metropolis algorithm. The optimiza- 
tion consists in the addition to the Metropolis al- 
gorithm of a second Markov chain in the temper- 
ature parameter (3 = J\/T. The new distribution 
law is in this case written as 



L = 24, T = 3.2, 5T = 0.2 
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FIG. 1: Overlap of energy distributions 
for adjacent replicas. 



P{{*}, {Pa}) ~ ex V (-p a H({a}) + g 



(5) 



g a is a certain constant. To each f3 a there corresponds its own g a . The probability of 
changing the dynamic parameter (3 a will obey an exp (-S), law, where 



S=((3' a -(3 a )H({a}) + (g' a -g a ). 



(6) 



In the case of the simultaneous simulation of the replicas of the system for each f3 a with 
a — . . . N, the change in the temperature can depend only on the nearest values of a; i.e., 
we have /3' a = /3 a ±i- The probability of a transition between the states that are determined 
by adjacent temperatures will be written as exp (—AS), where AS is assigned by the 
expression 

AS — S' — S — (E Pa - Ep a+1 )((3 a+1 - (3 a ). (7) 

According to this algorithm, the temperature of the system can in the process of simulation 
both decrease (annealing of the system) and increase, which makes it possible for the system 
to surmount high potential barriers. 

The basic criterion for the application of the parallel-tempering method is the overlap 
of the energy-distribution functions of the replicas of the system for adjacent temperatures 
(Fig. [T]). Since the energy variance approaches zero with decreasing temperature, the selec- 
tion of equidistant temperatures is not justified; the interval between adjacent temperatures 



ST = T a+ i — T a also should decrease. This dependence leads to additional difficulties upon 
the simulation of the behavior of these systems in the region of low temperatures, since 
the initial maximum temperature of simulation Tn should be sufficiently large to surmount 
the potential barriers which separate the local energy minima. Meanwhile, the number 
of replicas of the system (at different temperatures) considered in simulation must ensure 
the mechanism of temperature exchange. Under the conditions of a finite number of the 
simulated replicas, which is limited by the computational resources of the researcher, it is 
necessary to solve the problem of an optimum choice of the number of temperatures to 
be considered. Theoretically, the most optimum selection of temperatures, according to 
I5I . , is such at which the probability of a transition of the replica to a new temperature 
is constant for the entire set of the replicas: 



As an optimum sequence of temperatures, we can use the geometric progression T a+ i/T a = 
const with: 



where R = ^/T N /T . 

The choice of the temperature sequence in this form ensures the fulfillment of equality 
(JBJ). The optimum value of the probability of the replica transition to a new temperature 
upon the realization of the parallel-tempering algorithm was determined [17| to be P ~ 0.23. 

However, in contrast to the Ising spin glasses, in which the temperature of transition into 
the spin-glass state is Tf/J ph 1, in the dilute antiferromagnetic Ising model with the effects 
of random magnetic fields at a spin concentration p = 0.5 the temperature of the phase 
transition into the mixed state is T m pa 5 in the units of exchange integral J\. As a result, 
the application of the parallel-tempering method to the simulation of such a system requires 
the use of a considerably wider temperature range. Thus, to guarantee obtaining a stable 
initial equilibrium state from which it is better to perform simulation, the initial maximum 
temperature Tjy must be selected to be considerably higher than the temperature T mix ; only 
after this it is possible to investigate the phase transition into the mixed state near T m i x 
in the vicinity of which there already appears a set of metastable states with an extremely 
slow dynamics of establishing equilibrium. On the other hand, for studying the asymptotic 
approach to the ground spin-glass state of the system it is desirable that the minimum 




(8) 




(9) 



temperature T be as close as possible to T = 0. All this imposes new restrictions on the 
realization of the parallel-tempering algorithm for studying disordered antiferromagnets and 
for the selection of the sequence of temperatures for the replicas to be simulated. 

Note also that for the realization of the condition of the equal probability of the tem- 
perature exchange between the adjacent replicas it is necessary to have information on the 
temperature dependence of the energy of the system E(T); for this reason, the procedure 
of finding an optimum set of temperature points requires a preliminary simulation of the 
system. As the first approximation, a sequence of temperatures in the form of a geometric 
progression ([9]) can be chosen. The subsequent simulation of the system for the selected 
temperatures gives an idea about the E(T) dependence, though also only in the first ap- 
proximation. At the next step, a new set of temperatures f3' a is determined according to 
condition (jSJ) to yield 

(E{(5' a ) - E((3' a+1 )) (f3' a+1 - #,) ~ const . (10) 

The new set of temperatures obtained in this case specifies new (more optimum) pa- 
rameters for the simulation of the disordered system. In turn, the thus obtained replica- 
temperature distribution is also an approximation (second-order), which can be used for the 
subsequent iterative search for the optimum set of temperatures. 

The inconvenience of this approach is the need to, each time anew, conduct the complete 
simulation of the system at each step of the search for an optimum sequence of temperatures. 
Nevertheless, this approach makes it possible to most effectively use the parallel-tempering 
algorithm in conducting numerical experiments on studying the low-temperature properties 
of spin-glass states. 

IV. RESULTS AND DISCUSSION 

In this work, we consider the antiferromagnetic Ising model with a spin concentration 
p = 0.5 (corresponding to the region of strong disordering) and a magnetic field h = 2. The 
simulation was conducted for a wide set of values of the linear dimensions of the cubic lattice 
(L = 8,16,24,32, and 40). 

For determining the temperature of the phase transition in the system and the intervals of 
existence of different phase states, we studied the temperature behavior of the heat capacity 
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FIG. 2: Temperature dependence of the heat 
capacity C for systems with a spin con- 
centration p = 0.50 and linear sizes L = 
8,16,24,32, and 40. 



FIG. 3: Temperature dependence of the 
spin-glass order parameter q for the systems 
with a spin concentration p = 0.50 and linear 
sizes L = 8, 16, 24, 32, and 40. 



of the system (Fig. [5]) for the lattices of the above-indicated dimensions. It is evident 
that in the temperature interval of T = 4.56.5 there is observed an anomalous increase in 
the heat capacity, which indicates the existence of a phase transition in the system. The 
size-dependent changes in the behavior of the heat capacity indicate the suppression of 
the energy fluctuations in the system in comparison with the typical second-order phase 
transitions into the antiferromagnetic state, and the bend in the temperature dependence 
C(T) at T mix = 5.13 for the lattice with L max = 40 is typical of the phase transition in spin 
glasses [3]. 

For obtaining equilibrium magnetic characteristics of the system, the initial states were 
selected in the paramagnetic phase. This choice is caused by the fact that, because of 
the presence of metastable states near the transition temperature and in the entire low- 
temperature phase, there arises a problem of reaching equilibrium initial configurations. 
The initial states obtained were used in the realization of the parallel-tempering method. 
Within the framework of this method, the initial set of temperatures was selected according 
to formula ([9]). The temperature dependence of the energy E(T) obtained in the simulation 
was used for refining the temperatures f3 a according to the principle of the equal probability 
of transition between adjacent temperatures. To achieve the equilibrium state of the system 
at each temperature, 10 4 Monte Carlo steps were used for the relaxation, with the rejection 
of half initial configurations. 



FIG. 4: Temperature dependence of the 
staggered (M stg , upper curves) and total 
(M, lower curves) magnetizations for sys- 
tems with a spin concentration p = 0.50 and 
linear sizes L = 8, 16, 24, 32, and 40. 



FIG. 5: Temperature dependence of staggered 
(M stg , upper curves) and total (M, lower curves) 
magnetizations for systems with a spin concen- 
tration p = 0.90 and linear sizes L = 24, 32, and 
40. In the inset, the temperature range of phase 
transition is shown for M s t g . 



Figures [3H display the temperature dependences obtained for the staggered (M stg ) and 
total (M) magnetizations and spin-glass order parameter (g) for the lattices with linear di- 
mensions from L = 8 to L = 40 averaged over 100 different impurity configurations. It is 
seen from the figures that all the values measured demonstrate a noticeable dependence on 
the sizes of the system. The strongest dependence on L is characteristic of the staggered 
magnetization, which for the systems of small sizes determines the prevailing magnetic or- 
dering of antiferromagnetic nature; the spin-glass order parameter in this case appears as a 
secondary parameter of ordering and repeats the temperature dependence of the staggered 
magnetization. The strong decrease in M stg with increasing L at noticeably smaller changes 
in the spin-glass order parameter indicates the predominance of spin-glass ordering in the 
system with L > 24 and the appearance, in the low-temperature phase, of a mixed phase 
state consisting of antiferromagnetic and ferromagnetic domains surrounded by a spin-glass 
phase. Note that at a temperature T mix = 5.13 the temperature dependence of the total 
magnetization for the lattice with L = 40 reaches a maximum with a subsequent character- 
istic decrease in M(T) in the high-temperature region, as was observed in 131 ] . The results 
obtained indicate that, in the limit of L — > oo and T — > 0, a spin-glass ground state is 



realized in the system. 

To compare the observed features of the 
behavior of the magnetic characteristics of 
the model in the region of strong struc- 
tural disorder with the behavior of analo- 
gous characteristics in the region of weak 
dilution, we carried out an analogous study 
for the samples with a spin concentration 
p = 0.9. Figure [5] displays the temper- 
ature dependences of the staggered (M stg ) 
and total (M) magnetizations on the linear 
dimensions of the lattice L = 24, 32, and 
40. It is evident that no size dependence of 
these characteristics is observed in the low- 
temperature phase, except for the tempera- 
ture range close to the critical temperature 
T c = 8.73. The temperature behavior of 
M stg clearly indicates the existence of an an- 
tiferromagnetic ordering in the system with 
an the antiferromagnetic ground state under 
above conditions. 

To confirm the existence of a domain 
structure in the mixed phase state of the 
strongly disordered antiferromagnet with 
random fields, we realized, within the frame- 
work of the statistical parallel-tempering 
method, a study of the temperature depen- 
dence of the local values of magnetic char- 
acteristics for cubic blocks of size 6 = 5, 10, 
and 20 into which the lattice under consid- 
eration (with a maximum size L = 40) was 
divided. These dependences are shown in Fi, 
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FIG. 6: Temperature dependence of the local 
values of the staggered magnetization M stg for 
blocks with sizes 6 = 5, 10, and 20. 
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FIG. 7: Temperature dependence of the local 
values of the total magnetization M for blocks 
with sizes 6 = 5, 10, and 20. 




0,0-1 1 , 1 , 1 , 1 , 1 , 1 , 1 , r 

1 2 3 4 5 6 7 

T 



FIG. 8: Temperature dependence of the local 
values of the spin-glass order parameter q for 
blocks with sizes 6 = 5, 10, and 20. 
,s. HE! solid curves correspond to the depen- 



dence of the average values of M stg (T), M(T), and q(T) for the entire lattice with L = 40; 
curves 1 and 1', 2 and 2', and 3 and 3' correspond to minimum and maximum values of 
these values for the blocks with 6 = 5, 10, and 20, respectively. An analysis of these figures 
and the data on the entire totality of blocks shows that as the temperature decreases, the 
sizes of typical antiferromagnetic domains decrease from l a ~ 20 to la ~ 10, and the sizes 
of ferromagnetic domains, from If ~ 10 to If ~ 5 with an increase in the volume of the 
spin-glass phase, until at T = there is realized a spin-glass ground state. 

Conclusions 

Thus, by applying the numerical parallel-tempering method for studying the low- 
temperature behavior of the three-dimensional disordered antiferromagnetic model with 
random magnetic fields, it was clearly shown that in weakly disordered systems there is 
realized an antiferromagnetic ordered state, whereas in the region of strong structural dis- 
order the effects of random magnetic fields lead to the realization of a new phase state in 
the system, which is characterized by a complex domain structure consisting of antiferro- 
magnetic and ferromagnetic domains separated by regions of a spin-glass phase with the 
realization of a spin-glass ground state. 
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